The smt_im data

• 894 subjects in Cleveland-Elyria with bmi and no history of diabetes (missing values singly imputed: assume MAR)
• All subjects have hx_diabetes (all 0), and are located in the MMSA labeled Cleveland-Elyria.
• See Course Notes Chapter on BRFSS SMART data for variable details
• Appendix provides details on data development.

The Five Variables We’ll Use Today

9 variables in the data but we’ll use only these 5 today.

Data Load

smt_im <- read_rds("c03/data/smt_im.Rds") |>
  select(ID, bmi, exerany, health, fruit_day, everything())

smt_im

# A tibble: 894 × 9
   ID      bmi exerany health fruit_day inc_imp drinks_wk female race_eth       
   <chr> <dbl>   <dbl> <fct>      <dbl>   <dbl>     <dbl>  <dbl> <fct>          
 1 2      23.0       0 E           4      86865      0         1 White non-Hisp…
 2 3      26.9       1 G           3      20904      0         1 Other race non…
 3 4      26.5       1 G           2      46946      4.67      1 White non-Hisp…
 4 5      24.2       1 G           0.57   58311      0.93      0 White non-Hisp…
 5 7      23.0       1 G           2       2318      2         0 White non-Hisp…
 6 8      28.4       1 VG          1      79667      0         0 Other race non…
 7 9      30.1       1 F           0.23   47880      0         0 Black non-Hisp…
 8 10     19.8       1 E           0.77  100136      0.47      1 White non-Hisp…
 9 11     27.2       1 E           0.71   73145      0         0 White non-Hisp…
10 12     24.6       1 E           1.07   76917      0         1 Other race non…
# ℹ 884 more rows

Checking our Data

Are there any missing values?

smt_im |> n_miss() 

[1] 0

Does each row have a unique ID value?

identical(nrow(smt_im), n_distinct(smt_im$ID)) 

[1] TRUE

Range and Level Checks?

data_codebook(smt_im |> select(bmi, exerany, health, fruit_day)) 

select(smt_im, bmi, exerany, health, fruit_day) (894 rows and 4 variables, 4 shown)

ID | Name      | Type        | Missings |     Values |           N
---+-----------+-------------+----------+------------+------------
1  | bmi       | numeric     | 0 (0.0%) | [13.3, 63] |         894
---+-----------+-------------+----------+------------+------------
2  | exerany   | numeric     | 0 (0.0%) |          0 | 232 (26.0%)
   |           |             |          |          1 | 662 (74.0%)
---+-----------+-------------+----------+------------+------------
3  | health    | categorical | 0 (0.0%) |          E | 148 (16.6%)
   |           |             |          |         VG | 324 (36.2%)
   |           |             |          |          G | 274 (30.6%)
   |           |             |          |          F | 113 (12.6%)
   |           |             |          |          P |  35 ( 3.9%)
---+-----------+-------------+----------+------------+------------
4  | fruit_day | numeric     | 0 (0.0%) |    [0, 10] |         894
------------------------------------------------------------------

Our outcome, bmi

ggplot(smt_im, aes(x = bmi)) +
  geom_histogram(binwidth = 2, col = "azure", fill = "coral")

Key predictors: exerany, health

smt_im |> tabyl(exerany, health) |> 
  adorn_totals(where = c("row", "col")) |>
  gt() |> tab_options(table.font.size = 28)

Here, it doesn’t matter much whether we store the 1/0 in exerany as numeric or as a two-level factor in R. For binary variables, sometimes the numeric version will be more useful and sometimes a factor will be more useful.

Our covariate, fruit_day

We are mostly interested in whether accounting for the quantitative covariate fruit_day changes the modeled association of our key predictors with bmi.

• Sometimes we center such a covariate (subtracting its mean.)

smt_im <- smt_im |>
  mutate(fruit_c = fruit_day - mean(fruit_day))

• Why? So that we can easily plug in the covariate’s mean (which will now be 0) when making predictions.

Did we center fruit_day properly?

Here’s a little “sanity check”:

df_stats(~ fruit_day + fruit_c, data = smt_im) |> gt() |>
  fmt_number(columns = min:sd, decimals = 3) |>
  tab_options(table.font.size = 24)

The df_stats() function comes from the mosaic package, and can be used to apply favstats() to multiple variables at once.

Modeling Plan

1. Split smt_im into training and testing samples.
2. Predict bmi using exer_any and health
   • (fit1): without an interaction between the predictors
   • (fit2): and then with an interaction term
3. (fit3): Then add in our (centered) covariate, fruit_c.
4. (fit4): Use a quadratic polynomial in fruit_c.
5. Assess all four models in training and testing samples.

Splitting the Sample

We’ll partition our data set using some tools from the rsample package, into:

• a training sample containing 75% of the data
• a testing sample containing the remaining 25%

set.seed(432)    ## for future replication

smt_im_split <- initial_split(smt_im, prop = 3/4)

train_smt_im <- training(smt_im_split)
test_smt_im <- testing(smt_im_split)

c(nrow(smt_im), nrow(train_smt_im), nrow(test_smt_im))

[1] 894 670 224

Modeling Plan

• Predict bmi using exer_any and health
  • (fit1): without an interaction between the predictors
  • (fit2): and then with an interaction term
• (fit3): Then add in our (centered) covariate, fruit_c.
• (fit4): Use a quadratic polynomial in fruit_c.

Tukey’s ladder of power transformations

• Used in combination with a Box-Cox plot from car
• Requires the \(y\) variable to be strictly positive.
  • If desirable, you can add any constant to \(y\) or multiply \(y\) by any constant, before or after the transformation.
• Be sure you can back out of the transformation:
• \(exp(log(y)) = y\), \(\sqrt{y^2} = y\), \(1/\frac{1}{y} = y\), \([1 / (\frac{1}{\sqrt{y}})]^2 = y\)

Consider transforming bmi?

m0 <- lm(bmi ~ exerany + health + fruit_c, data = train_smt_im)
boxCox(m0)

Should we transform bmi? (n = 670)

p1a <- ggplot(train_smt_im, aes(x = bmi)) + 
    geom_histogram(col = "navy", fill = "cyan", bins = 20)

p1b <- ggplot(train_smt_im, aes(x = bmi, y = "bmi")) + 
    geom_boxplot(col = "navy", fill = "cyan") + labs(y = "") +
    stat_summary(geom = "point", col = "red", fun = mean, size = 3)

p2a <- ggplot(train_smt_im, aes(x = 1/sqrt(bmi))) + 
    geom_histogram(col = "navy", fill = "gold", bins = 20)

p2b <- ggplot(train_smt_im, aes(x = 1/sqrt(bmi), y = "1/sqrt(bmi)")) +
    geom_boxplot(col = "navy", fill = "gold") + labs(y = "") +
    stat_summary(geom = "point", col = "red", fun = mean, size = 3)

p3a <- ggplot(train_smt_im, aes(x = 1/bmi)) + 
    geom_histogram(col = "navy", fill = "green", bins = 20)

p3b <- ggplot(train_smt_im, aes(x = 1/bmi, y = "1/bmi")) + 
    geom_boxplot(col = "navy", fill = "green") + labs(y = "") +
    stat_summary(geom = "point", col = "red", fun = mean, size = 3)

(p1a + p1b) / (p2a + p2b) / (p3a + p3b)

Re-scaling the transformation

To ease interpretation of coefficients, I sometimes scale an outcome transformation so that its values fall in (10, 100), rather than between 0 and 1.

bind_rows( favstats(~ 1/sqrt(bmi), data = train_smt_im),
           favstats(~ 100/sqrt(bmi), data = train_smt_im)) |>
  mutate(outcome = c("1/sqrt(bmi)", "100/sqrt(bmi)")) |> 
  relocate(outcome) |>
  gt() |> fmt_number(columns = min:sd, decimals = 3) |> 
  tab_options(table.font.size = 24)

Shape doesn’t change

Means by exerany and health

summaries_1 <- train_smt_im |>
    group_by(exerany, health) |>
    summarise(n = n(), mean = mean(100/sqrt(bmi)), 
              stdev = sd(100/sqrt(bmi)))
summaries_1 

# A tibble: 10 × 5
# Groups:   exerany [2]
   exerany health     n  mean stdev
     <dbl> <fct>  <int> <dbl> <dbl>
 1       0 E         18  19.2  1.22
 2       0 VG        55  19.5  1.90
 3       0 G         60  18.6  2.22
 4       0 F         32  17.5  2.43
 5       0 P          9  17.4  2.25
 6       1 E         92  19.9  1.64
 7       1 VG       189  19.6  1.74
 8       1 G        150  18.8  1.91
 9       1 F         49  19.5  1.94
10       1 P         16  19.3  2.79

Code for Interaction Plot

ggplot(summaries_1, aes(x = health, y = mean, col = factor(exerany))) +
  geom_line(aes(group = factor(exerany)), linewidth = 2) +
  scale_color_viridis_d(option = "C", end = 0.5) +
  labs(title = "Observed Means of 100/sqrt(BMI)",
       subtitle = "by Exercise and Overall Health")

• Note the use of factor here since the exerany variable is in fact numeric, although it only takes the values 1 and 0.
  • Sometimes it’s helpful to treat 1/0 as a factor, and sometimes not.
• Where is the evidence of serious non-parallelism (if any) in the plot (see next slide) that results from this code?

Create our transformed outcome

We’ll want to do this in both our training and test samples.

train_smt_im <- train_smt_im |> mutate(bmi_tr = 100 / sqrt(bmi))

test_smt_im <- test_smt_im |> mutate(bmi_tr = 100 / sqrt(bmi))

Model fit1 without interaction

fit1 <- lm(bmi_tr ~ exerany + health, data = train_smt_im)

Using the tidy() function from broom:

tidy(fit1, conf.int = TRUE, conf.level = 0.90) |> 
  gt() |> fmt_number(columns = estimate:conf.high, decimals = 3) |>
  tab_options(table.font.size = 24)

Model Parameters for fit1

model_parameters(fit1, ci = 0.90) 

Parameter   | Coefficient |   SE |         90% CI | t(664) |      p
-------------------------------------------------------------------
(Intercept) |       19.31 | 0.23 | [18.92, 19.69] |  82.84 | < .001
exerany     |        0.57 | 0.17 | [ 0.29,  0.85] |   3.33 | < .001
health [VG] |       -0.20 | 0.22 | [-0.56,  0.16] |  -0.90 | 0.367 
health [G]  |       -0.96 | 0.23 | [-1.34, -0.59] |  -4.23 | < .001
health [F]  |       -0.98 | 0.28 | [-1.45, -0.51] |  -3.46 | < .001
health [P]  |       -1.08 | 0.43 | [-1.79, -0.38] |  -2.53 | 0.012 

Model Parameters for fit1 (with gt())

Reformatting with gt()…

model_parameters(fit1, ci = 0.90) |> 
  gt() |> fmt_number(columns = -c(CI, df_error), decimals = 3) |>
  tab_options(table.font.size = 24)

• Effect of exerany on \(100/\sqrt{bmi}\)?
• Effect of health = Fair instead of Excellent?

How well does fit1 fit the training data?

n_obs(fit1)

[1] 670

model_performance(fit1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
2786.766 | 2786.935 | 2818.317 | 0.069 |     0.062 | 1.916 | 1.925

glance(fit1) |> 
    select(r.squared, adj.r.squared, sigma, nobs, 
           df, df.residual, AIC, BIC) |> 
  gt() |> fmt_number(columns = r.squared:sigma, decimals = 3) |>
  fmt_number(columns = AIC:BIC, decimals = 1) |>
  tab_options(table.font.size = 24)

Tidied ANOVA for fit1

tidy(anova(fit1)) |> gt() |> 
  fmt_number(columns = sumsq:statistic, decimals = 2) |>
  fmt_number(columns = p.value, decimals = 4) |>
  tab_options(table.font.size = 24)

Model Checks

We’ll be checking assumptions related to:

• linearity
• homoscedasticity (constant variance)
• influential observations (outliers, leverage and influence)
• whether the residuals follow a Normal distribution
• collinearity (variance inflation factor)
• and a posterior predictive check of our predictions

My slides and check_model()

When building a regular HTML file, I would just use:

check_model(fit1, detrend = FALSE)

with #| fig-height: 9 at the start of the code chunk so that the plots are taller than the default height (thus easier to read) but I will split out the plots for slides.

Checking model fit1 (n = 670)

check_model(fit1, check = c("linearity", "homogeneity"))

Checking model fit1 (n = 670)

check_model(fit1, check = c("outliers", "qq"), detrend = FALSE)

Checking model fit1 (n = 670)

check_model(fit1, check = c("pp_check", "vif"))

Adding the interaction term to fit1

fit2 <- lm(bmi_tr ~ exerany * health, data = train_smt_im)

• How do our models compare on fit to the training data?

bind_rows(glance(fit1), glance(fit2)) |>
  mutate(mod = c("fit1", "fit2")) |>
  select(mod, r.sq = r.squared, adj.r.sq = adj.r.squared, 
         sigma, nobs, df, df.res = df.residual, AIC, BIC) |> 
  gt() |> fmt_number(columns = r.sq:sigma, decimals = 3) |>
  fmt_number(columns = AIC:BIC, decimals = 1) |>
  tab_options(table.font.size = 24)

ANOVA for the fit2 model

tidy(anova(fit2)) |> gt() |> 
  fmt_number(columns = sumsq:statistic, decimals = 2) |>
  fmt_number(columns = p.value, decimals = 4) |>
  tab_options(table.font.size = 20)

ANOVA test comparing fit1 to fit2

anova(fit1, fit2)

Analysis of Variance Table

Model 1: bmi_tr ~ exerany + health
Model 2: bmi_tr ~ exerany * health
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
1    664 2459.8                                
2    660 2392.3  4    67.544 4.6586 0.001029 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit2 coefficients

tidy(fit2, conf.int = TRUE, conf.level = 0.90) |>
  gt() |> fmt_number(columns = estimate:conf.high, decimals = 3) |>
  tab_options(table.font.size = 20)

Interpreting the fit2 model

• How do we interpret effect sizes here? It depends…

Interpreting the fit2 model

• Effect of exerany on predicted \(100/\sqrt{bmi}\)?
  • If health = Excellent, effect is +0.72
  • If health = Fair, effect is (0.72 + 1.29) = +2.01
• Effect of health = Fair instead of Excellent?
  • If exerany = 0 (no), effect is -1.73
  • If exerany = 1 (yes), effect is (-1.73 + 1.29) = -0.44

Checking model fit2 (n = 670)

check_model(fit2, check = c("linearity", "homogeneity"))

Checking model fit2 (n = 670)

check_model(fit2, check = c("outliers", "qq"), detrend = FALSE)

Checking model fit2 (n = 670)

check_model(fit2, check = c("pp_check", "vif"))

Add fruit_c to fit1

fit3 <- lm(bmi_tr ~ fruit_c + exerany + health, data = train_smt_im)

• How well does this model fit the training data?

bind_rows(glance(fit1), glance(fit3)) |>
  mutate(mod = c("fit1", "fit3")) |>
  select(mod, r.sq = r.squared, adj.r.sq = adj.r.squared, 
         sigma, df, df.res = df.residual, AIC, BIC) |> 
  gt() |> fmt_number(columns = r.sq:sigma, decimals = 3) |>
  fmt_number(columns = AIC:BIC, decimals = 1) |>
  tab_options(table.font.size = 24)

ANOVA for the fit3 model

tidy(anova(fit3)) |> gt() |> 
  fmt_number(columns = sumsq:statistic, decimals = 2) |>
  fmt_number(columns = p.value, decimals = 4) |>
  tab_options(table.font.size = 24)

fit3 coefficients

tidy(fit3, conf.int = TRUE, conf.level = 0.90) |>
  gt() |> fmt_number(columns = estimate:conf.high, decimals = 3) |>
  tab_options(table.font.size = 24)

Checking model fit3 (n = 670)

check_model(fit3, detrend = FALSE, check = c("pp_check", "qq"))

Include the interaction term?

fit4 <- lm(bmi_tr ~ fruit_c + exerany * health, 
          data = train_smt_im)

ANOVA for the fit4 model

tidy(anova(fit4)) |> gt() |> 
  fmt_number(columns = sumsq:statistic, decimals = 2) |>
  fmt_number(columns = p.value, decimals = 4) |>
  tab_options(table.font.size = 20)

fit4 coefficients

tidy(fit4, conf.int = TRUE, conf.level = 0.90) |>
  gt() |> fmt_number(columns = estimate:conf.high, decimals = 3) |>
  tab_options(table.font.size = 18)

ANOVA: Compare fit3 & fit4

anova(fit3, fit4)

Analysis of Variance Table

Model 1: bmi_tr ~ fruit_c + exerany + health
Model 2: bmi_tr ~ fruit_c + exerany * health
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)    
1    663 2437.2                                 
2    659 2365.0  4    72.145 5.0257 0.000539 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Checking model fit4 (n = 670)

check_model(fit4, detrend = FALSE, check = c("pp_check", "qq"))

Which of the four models fits best?

In the training sample, our model results (note ordering):

• Adjusted \(R^2\), \(\sigma\) and AIC all improve as we move down this table. BIC likes fit1 and fit3.
• The training sample is the data our models have already seen, so we should be cautious.

Comparison Plot: In-Sample Performance

plot(compare_performance(fit1, fit2, fit3, fit4))

What does augment() give us?

fit1_test_aug <- augment(fit1, newdata = test_smt_im) 
fit1_test_aug |> select(ID, bmi_tr, bmi, .fitted, .resid, health, exerany) |>
  slice(198:202) |> gt() |> 
  fmt_number(columns = bmi_tr:.resid, decimals = 2) |>
  tab_options(table.font.size = 20)

Here, .fitted = predicted bmi_tr and .resid = bmi_tr - .fitted.

Back-Transformation of bmi_tr

Our models predict bmi_tr = \(100/\sqrt{bmi}\), but we want to predict bmi. How do we convert predicted \(100/\sqrt{bmi}\) to predicted bmi?

\[ 1 / \left(\frac{100}{\sqrt{bmi}}\right) = \sqrt{bmi} / 100, \\ \mbox{so } 100 / \left(\frac{100}{\sqrt{bmi}}\right) = \sqrt{bmi}, \\ \mbox{and so } \left[100 / \left(\frac{100}{\sqrt{bmi}}\right)\right]^2 = bmi \]

augment() with results for bmi

We use \((\frac{100}{.fitted})^2\) for predicted bmi, then errors are bmi_res = observed bmi - predicted bmi.

fit1_test_aug <- augment(fit1, newdata = test_smt_im) |> 
  mutate(bmi_fit = (100/.fitted)^2, bmi_res = bmi - bmi_fit)

fit1_test_aug |> select(ID, bmi, bmi_fit, bmi_res, 
      bmi_tr, .fitted, .resid, exerany, health, fruit_c) |>
  slice(5:6) |> gt() |> 
  fmt_number(columns = c(bmi:.resid, fruit_c), decimals = 2) |>
  tab_options(table.font.size = 24)

Augment all four models so far…

fit1_test_aug <- augment(fit1, newdata = test_smt_im) |>
  mutate(bmi_fit = (100/.fitted)^2, bmi_res = bmi - bmi_fit)

fit2_test_aug <- augment(fit2, newdata = test_smt_im) |>
  mutate(bmi_fit = (100/.fitted)^2, bmi_res = bmi - bmi_fit)

fit3_test_aug <- augment(fit3, newdata = test_smt_im) |>
  mutate(bmi_fit = (100/.fitted)^2, bmi_res = bmi - bmi_fit)

fit4_test_aug <- augment(fit4, newdata = test_smt_im) |>
  mutate(bmi_fit = (100/.fitted)^2, bmi_res = bmi - bmi_fit)

Four Key Error Summaries

We’ll look at all four of these summaries when we do linear regression, usually.

• Mean absolute prediction error (MAPE)
• Maximum absolute prediction error (Max. Error)
• Square root of mean squared prediction error (RMSPE)
• Squared correlation of observed and predicted bmi (validated \(R^2\))

Key Summaries for fit1

fit1_esum <- fit1_test_aug |>
  summarise(MAPE = mean(abs(bmi_res)),
            Max_E = max(abs(bmi_res)),
            RMSPE = sqrt(mean(bmi_res^2)),
            Val_R2 = cor(bmi, bmi_fit)^2) |>
  mutate(Model = "fit1")

fit1_esum

# A tibble: 1 × 5
   MAPE Max_E RMSPE Val_R2 Model
  <dbl> <dbl> <dbl>  <dbl> <chr>
1  4.31  21.1  5.63 0.0752 fit1 

• I built the key summaries for fit2, fit3 and fit4 in the same way (included in code, not shown in slides.)

Compare Models in Test Sample

bind_rows(fit1_esum, fit2_esum, fit3_esum, fit4_esum) |>
  relocate(Model) |> gt() |> fmt_number(decimals = 3) |>
  tab_options(table.font.size = 24)

Next up…

Basics of logistic regression fitting and evaluation

• What if we have a binary (yes/no or 1/0) outcome?
• Predict “whether or not BMI < 30”, rather than BMI?
  • A linear probability model as a first idea
  • Using glm() rather than lm() to get a logistic model
  • Coefficients as log(odds ratios)
  • Changes in how we measure the model’s performance
  • Changes in the assumptions we make

Creating Today’s Data Set

url1 <- "https://raw.githubusercontent.com/THOMASELOVE/432-data/master/data/smart_ohio.csv"

smart_ohio <- read_csv(url1)

smt <- smart_ohio |>
    filter(hx_diabetes == 0, mmsa == "Cleveland-Elyria",
           complete.cases(bmi)) |>
    select(bmi, inc_imp, fruit_day, drinks_wk, 
           female, exerany, genhealth, race_eth, 
           hx_diabetes, mmsa, SEQNO) |>            
    mutate(across(where(is.character), as_factor)) |>
    mutate(ID = as.character(SEQNO - 2017000000)) |>
    relocate(ID)

Codebook for useful smt variables (1)

• 894 subjects in Cleveland-Elyria with bmi and no history of diabetes

Codebook for useful smt variables (2)

• 894 subjects in Cleveland-Elyria without diabetes

• plus ID, SEQNO, hx_diabetes (all 0), MMSA
• See Course Notes Chapter 6 on BRFSS SMART data

Basic Data Summaries

Available approaches include:

• data_codebook() from datawizard in easystats
• Hmisc package’s describe(), or
• summary()

all of which can work nicely in an HTML presentation, but none of them fit well on a slide.

Histogram of each quantity

Warning: Removed 120 rows containing non-finite values

p1 <- ggplot(smt, aes(x = bmi)) + 
    geom_histogram(fill = "navy", col = "white", bins = 20)
p2 <- ggplot(smt, aes(x = inc_imp)) + 
    geom_histogram(fill = "forestgreen", col = "white", bins = 20)
p3 <- ggplot(smt, aes(x = fruit_day)) + 
    geom_histogram(fill = "tomato", col = "white", bins = 20)
p4 <- ggplot(smt, aes(x = drinks_wk)) + 
    geom_histogram(fill = "royalblue", col = "white", bins = 20)

(p1 + p2) / (p3 + p4)

Binary variables in raw smt

smt |> tabyl(female, exerany) |> adorn_title()

        exerany        
 female       0   1 NA_
      0      95 268  20
      1     128 361  22

• female is based on biological sex (1 = female, 0 = male)
• exerany comes from a response to “During the past month, other than your regular job, did you participate in any physical activities or exercises such as running, calisthenics, golf, gardening, or walking for exercise?” (1 = yes, 0 = no, don’t know and refused = missing)
• Any signs of trouble here?

Multicategorical genhealth in raw smt

smt |> tabyl(genhealth)

   genhealth   n     percent valid_percent
 1_Excellent 148 0.165548098    0.16573348
      3_Good 274 0.306487696    0.30683091
  2_VeryGood 324 0.362416107    0.36282195
      4_Fair 112 0.125279642    0.12541993
      5_Poor  35 0.039149888    0.03919373
        <NA>   1 0.001118568            NA

• The variable is based on “Would you say that in general your health is …” using the five specified categories (Excellent -> Poor), numbered for convenience after data collection.
• Don’t know / not sure / refused treated as missing.
• How might we manage this variable?

Changing the levels for genhealth

smt <- smt |>
    mutate(health = 
               fct_recode(genhealth,
                          E = "1_Excellent",
                          VG = "2_VeryGood",
                          G = "3_Good",
                          F = "4_Fair",
                          P = "5_Poor"),
           health = fct_relevel(health, "E", "VG", "G", "F", "P"))

Might want to run a sanity check here, just to be sure…

Checking health vs. genhealth

smt |> tabyl(genhealth, health) |> adorn_title()

             health                   
   genhealth      E  VG   G   F  P NA_
 1_Excellent    148   0   0   0  0   0
      3_Good      0   0 274   0  0   0
  2_VeryGood      0 324   0   0  0   0
      4_Fair      0   0   0 112  0   0
      5_Poor      0   0   0   0 35   0
        <NA>      0   0   0   0  0   1

• OK. We’ve adjusted the order to something more sensible, retained the missing value, and we have much shorter labels.

Multicategorical race_eth in raw smt

smt |> count(race_eth)

# A tibble: 6 × 2
  race_eth                     n
  <fct>                    <int>
1 White non-Hispanic         646
2 Other race non-Hispanic     22
3 Black non-Hispanic         167
4 Multiracial non-Hispanic    19
5 Hispanic                    27
6 <NA>                        13

“Don’t know”, “Not sure”, and “Refused” were treated as missing.

• What is this variable actually about?
• What is the most common thing people do here?

What is the question you are asking?

Collapsing race_eth levels might be rational for some questions.

• We have lots of data from two categories, but only two.
• Systemic racism affects people of color in different ways across these categories, but also within them.

Is combining race and Hispanic/Latinx ethnicity helpful?

It’s hard to see the justice in collecting this information and not using it in as granular a form as possible, though this leaves some small sample sizes. There is no magic number for “too small a sample size.”

• Most people identified themselves in one category.
• These data are not ordered, and (I’d argue) ordering them isn’t helpful.
• Regression models are easier to interpret, though, if the “baseline” category is a common one.

Resorting the factor for race_eth

Let’s sort all five levels, from most observations to least…

smt <- smt |>
    mutate(race_eth = fct_infreq(race_eth))

smt |> tabyl(race_eth)

                 race_eth   n    percent valid_percent
       White non-Hispanic 646 0.72259508    0.73325766
       Black non-Hispanic 167 0.18680089    0.18955732
                 Hispanic  27 0.03020134    0.03064699
  Other race non-Hispanic  22 0.02460850    0.02497162
 Multiracial non-Hispanic  19 0.02125280    0.02156640
                     <NA>  13 0.01454139            NA

• Not a perfect solution, certainly, but we’ll try it out.

“Cleaned” Data and Missing Values

smt <- smt |>
    select(ID, bmi, inc_imp, fruit_day, drinks_wk, 
           female, exerany, health, race_eth)

miss_var_summary(smt)

# A tibble: 9 × 3
  variable  n_miss pct_miss
  <chr>      <int>    <num>
1 inc_imp      120   13.4  
2 exerany       42    4.70 
3 fruit_day     41    4.59 
4 drinks_wk     39    4.36 
5 race_eth      13    1.45 
6 health         1    0.112
7 ID             0    0    
8 bmi            0    0    
9 female         0    0    

Single Imputation with mice

smt_im <- mice(smt, m = 1, seed = 20250121, print = FALSE) |>
  complete() |>
  tibble()

prop_miss_case(smt_im)

[1] 0

dim(smt_im)

[1] 894   9

Saving the tidied data

Let’s save both the unimputed and the imputed tidy data as R data sets.

write_rds(smt, "c03/data/smt.Rds")

write_rds(smt_im, "c03/data/smt_im.Rds")

To reload these files, we’ll use read_rds().

• The main advantage here is that we’ve saved the whole R object, including all characteristics that we’ve added since the original download.